#### Replication for Tables and Figures in Paper ####
# Dylan Potts, EUI

#Packages
library(pacman)
p_load(tidyverse, estimatr, stargazer, fixest, fwildclusterboot, jtools, dotwhisker, marginaleffects,
       sjPlot, sjmisc, lfe, modelsummary, car, KalbeExtra, broom.mixed)

# Load Datasets- these can also be found on dataverse
load("Data/datasets/msr_army_trim.RData")
# Case Study datasets
# first the general dataset
load("Data/datasets/monthly_case.RData")
# Next the monthly_dataset
load("Data/datasets/types_case.RData")
# Donations dataset
load("Data/datasets/donations.RData")


# Table 1

hm1 <- felm(rh_cm ~ irish*fourties | 0 | 0 | comp_id, data = (msr_army_trim %>%  filter(american ==0)))
hm2 <- felm(rh_cm ~ irish*fourties | comp_id | 0 | comp_id, data = (msr_army_trim %>%  filter(american ==0)))
hm3 <- felm(rh_cm ~ irish*fourties | comp_id + enlisted_year | 0 | comp_id, data = (msr_army_trim %>%  filter(american ==0)))

stargazer::stargazer(hm1, hm2, hm3, 
                     keep = c("irish", "fourties", "irish:fourties"), 
                     type="latex",
                     title="Height Predicted by Cohorts Among Immigrant Soldiers",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "", "Y","Y" ), c("Enlistment Year FEs","", "","Y")),
                     covariate.labels = c("Irish", "1840s", "Irish x 1840s"),
                     dep.var.caption  = "Height at Enlistment", 
                     dep.var.labels   = "", 
                     notes.label = "Standard errors clustered by company",
                     column.separate = c(1,1),
                     header=FALSE, 
                     font.size = "small",
                     label = "tab:htest")

# Table 2

m1 <- felm(desert ~ irish*fourties | 0 | 0 | comp_id, data = msr_army_trim)
m2 <- felm(desert ~ irish*fourties | comp_id | 0 | comp_id, data = msr_army_trim)
m3 <- felm(desert ~ irish*fourties | comp_id + enlisted_year | 0 | comp_id, data = msr_army_trim)

stargazer::stargazer(m1, m2, m3, 
                     keep = c("irish", "fourties", "irish:fourties"), 
                     type="latex",
                     title="Cohorts and Desertion",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "", "Y","Y" ), c("Enlistment Year FEs","", "","Y")),
                     covariate.labels = c("Irish", "1840s", "Irish x 1840s"),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels   = "", 
                     notes.label = "Standard errors clustered by company",
                     column.separate = c(1,1),
                     header=FALSE,
                     font.size = "small",
                     label = "tab:main")


# Figure 1

famir3 <- feols(desert ~ younger  | famid + enlisted_year, data= (msr_army_trim%>% filter(irish==1) %>% 
                                                                    group_by(famid) %>% filter(n()>1) %>%
                                                                    filter(n_distinct(younger) ==2) %>% ungroup()))
famge3 <- feols(desert ~ youngerge | famid + enlisted_year, data= (msr_army_trim %>% filter(german==1) %>% group_by(famid) %>% filter(n()>1) %>%
                                                                     filter(n_distinct(younger) ==2) %>% ungroup()))
famam3 <- feols(desert ~ youngeram | famid + enlisted_year, data= (msr_army_trim %>% filter(american==1) %>% group_by(famid) %>% filter(n()>1) %>%
                                                                     filter(n_distinct(younger) ==2) %>% ungroup()))

sib_coef <- plot_coefs(famir3,famge3, famam3,
                       coefs = c("American"= "youngeram", 
                                 "German" = "youngerge",
                                 "Irish" = "younger"),
                       scale = TRUE, 
                       colors = c( "Aquamarine3", "Aquamarine3", "Aquamarine3" ), 
                       point.shape=FALSE, 
                       point.size=4
                       #groups = list(pane_1 = c("younger"),pane_2 = c("youngerge"), pane_3 = c("youngeram")))
)

sib_coef <- sib_coef + theme(legend.position="none")+coord_flip()
sib_coef <- sib_coef +
  ggtitle("")+
  xlab("Desertion")+
  ylab("Ethnic Group") +
  theme(axis.text.x = element_text(size = 20))


# Table 3

ig <- msr_army_trim %>%  filter(sibling==1) %>% 
  filter(american == 0) %>% 
  filter(english == 0) %>% 
  filter(!is.na(younger)) %>% 
  filter(freq==2) %>% 
  group_by(famid) %>%
  mutate(ff = length(famid)) %>% 
  filter(ff>=2) %>%
  mutate(y2 = sum(younger)) %>% 
  filter(y2==1) %>% 
  mutate(diff = ybirth[younger == 1] - ybirth) %>% 
  mutate(ad = mean(diff)*2) %>% 
  ungroup()
a <- msr_army_trim %>%  filter(sibling==1) %>% 
  filter(american == 1) %>% 
  filter(!is.na(younger)) %>% 
  filter(freq==2) %>% 
  group_by(famid) %>%
  mutate(ff = length(famid)) %>% 
  filter(ff>=2) %>%
  mutate(y2 = sum(younger)) %>% 
  filter(y2==1) %>% 
  mutate(diff = ybirth[younger == 1] - ybirth) %>% 
  mutate(ad = mean(diff)*2) %>% 
  ungroup()

aig <- ig %>%  bind_rows(a)

ilm <- felm(desert ~  younger * ad| 0 | 0 | comp_id, data= (ig %>% filter(irish==1)))
glm <- felm(desert ~  younger * ad| 0 | 0 | comp_id, data= (ig %>% filter(irish==0)))

stargazer::stargazer(ilm, glm,
                     keep = c("younger", "ad", "younger:ad"), 
                     type="latex",
                     title="Interaction Between Birth Order and Age Difference",
                     style="default",
                     keep.stat=c("n"),
                     covariate.labels = c("Younger", "Age Difference", "Younger x Age Difference"),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels   = "", 
                     column.labels = c("Irish Soldiers", "German Soldiers"),
                     notes.label = "",
                     column.separate = c(1,1),
                     header=FALSE, 
                     label = "tab:agediff")


# Table 4

mh1 <- felm(desert~ rh_cm*fourties | 0 | 0 | comp_id, data= (msr_army_trim %>% filter (irish ==1)))
mh2 <- felm(desert~ rh_cm*fourties | comp_id + enlisted_year | 0 | comp_id, data= (msr_army_trim %>% filter (irish ==1)))
mh3 <- felm(desert~ rh_cm*fourties | 0 | 0 | comp_id, data= (msr_army_trim %>% filter (irish ==0)))
mh4 <- felm(desert~ rh_cm*fourties | comp_id + enlisted_year | 0 | comp_id, data= (msr_army_trim %>% filter (irish ==0)))

stargazer::stargazer(mh1, mh2, mh3, mh4,
                     keep = c("rh_cm", "fourties", "rh_cm:fourties"), 
                     type="latex",
                     title="Interaction of Height and Being Born in the 1840s",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "", "Y","", "Y" ), c("Enlistment Year FEs","", "Y", "","Y")),
                     covariate.labels = c("Height", "1840s", "Height x 1840s"),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels = "",
                     column.labels   = c("Irish Soldiers", "Irish Soldiers", "Other Soldiers", "Other Soldiers"), 
                     notes.label = "Standard errors clustered by company",
                     header=FALSE, 
                     font.size = "small",
                     label = "tab:height")

# Table 5

msr_army_irish <- msr_army_trim %>% filter(irish==1) %>% mutate(works = log(1+works))



m_pop0 <- felm(desert~ exposureneg + ybirth + manual| enlisted_year + comp_id| 0 | comp_id, data=msr_army_irish)
m_pop1 <- felm(desert~ exposureneg + ybirth + manual| enlisted_year + comp_id + province| 0 | comp_id, data=msr_army_irish)
m_wor0 <- felm(desert~ works + ybirth + manual| enlisted_year + comp_id| 0 | comp_id, data=msr_army_irish)
m_wor1 <- felm(desert~ works + ybirth + manual| enlisted_year + comp_id + province| 0 | comp_id, data=msr_army_irish)
m_ind0 <- felm(desert~ findex + ybirth + manual| enlisted_year + comp_id| 0 | comp_id, data=msr_army_irish)
m_ind1 <- felm(desert~ findex + ybirth + manual| enlisted_year + comp_id + province| 0 | comp_id, data=msr_army_irish)


stargazer::stargazer(m_pop0, m_pop1, m_wor0, m_wor1, m_ind0, m_ind1, 
                     keep = c("exposureneg", "works", "findex"), 
                     type="latex",
                     title="Famine Intensity and Desertion",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "Y", "Y","Y", "Y","Y", "Y" ), 
                                    c("Enlistment Year FEs","Y", "Y", "Y","Y","Y", "Y"),
                                    c("Province FEs","", "Y", "","Y","", "Y")),
                     covariate.labels = c("Population Loss", "Log Public Works", "Famine Index"),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels = "",
                     column.labels   = "", 
                     notes.label = "Standard errors clustered by company",
                     header=FALSE, 
                     font.size = "small",
                     label = "tab:origins")
                     
# Table 6

mb <- msr_army_trim %>% filter(battles!=0) %>% 
  felm(desert~ irish*fourties | enlisted_year + comp_id | 0 | comp_id, data=.,)
mnb <- msr_army_trim %>% filter(battles==0) %>%
  felm(desert~ irish*fourties | enlisted_year + comp_id| 0 | comp_id, data=.,)

stargazer::stargazer(mb, mnb,
                     keep = c("irish", "fourties", "irish:fourties"), 
                     type="latex",
                     title="Cohorts With and Without Battles",
                     style="default",
                     keep.stat=c("n"),
                     covariate.labels = c("Irish", "1840s", "Irish x 1840s"),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels   = "", 
                     column.labels = c("Battles", "No Battles"),
                     notes.label = "Standard errors clustered by company",
                     column.separate = c(1,1),
                     header=FALSE, 
                     label = "tab:battles")

# Figure 3

df_date <- df_date %>% mutate(Nationality = ifelse(irish ==0, "German", "Irish"))

case_timeline <- qplot(x=month, y=desert, color=Nationality,
                       data=df_date, na.rm=TRUE)+
  labs(x="Time in the War", y="Number of Deserters", fill="Irish Troops")+
  theme_bw()
                     
# Figure 4

case_study <- case_study %>% mutate(Nationality = ifelse(irish ==0, "German", "Irish"))

case_bar<- ggplot(data=subset(case_study, !is.na(desert_type)), aes(x = factor(desert_type), y = desert, fill=Nationality)) + 
  geom_bar(stat='identity', position = "dodge")+
  labs(x="Type of Desertion", y="Number of Deserters", fill="Nationality")+
  theme_classic() +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14), 
        legend.text=element_text(size=14))


# Figure 5

ir_comp <- lm(desert~ enlisted_year + fourties*comp_ir, data=msr_army_irish)
ir_comp_mod <- plot_cap(ir_comp, condition = c("comp_ir", "fourties"))
ir_comp_mod <- ir_comp_mod + xlab("Share of a Company Which is Irish") +
  ylab("Predicted Desertion Rate") +
  labs(color = "Born in the 1840s")


# Table 7

manual1<- lm(manual ~ irish * fourties, data=msr_army_trim)

stargazer(manual1,
          keep = c("irish", "fourties", "irish:fourties"), 
          type="latex",
          title="The Role of Socioeconomic Status",
          style="default",
          keep.stat=c("n"),
          covariate.labels = c("Irish", "1840s", "Irish x 1840s"),
          dep.var.caption  = "", 
          dep.var.labels   = "",
          column.labels   = "Manual Worker" ,
          header=FALSE, 
          label="tab:manual")
            
                     
